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Abstract 

The susceptible-infectious-recovered (SIR) model describes the evolution of three species of in- 
dividuals which are subject to an infection and recovery mechanism. A susceptible S can become 
infectious with an infection rate /3 by an infectious /- type provided that both are in contact. The 
/- type may recover with a rate 7 and from then on stay immune. Due to the coupling between 
the different individuals, the model is nonlinear and out of equilibrium. We adopt a stochastic 
individual-based description where individuals are represented by nodes of a graph and contact 
is defined by the links of the graph. Mapping the underlying Master equation into a quantum 
formulation in terms of spin operators, the hierarchy of evolution equations can be solved exactly 
for arbitrary initial conditions on a linear chain. In case of uncorrelated random initial conditions 
the exact time evolution for all three individuals of the SIR model is given analytically. Depending 
on the initial conditions and reaction rates (5 and 7, the /-population may increase initially before 
decaying to zero. Due to fluctuations, isolated regions of susceptible individuals evolve and unlike 
in the standard mean-field SIR model one observes a finite stationary distribution of the S-type 
even for large population size. The exact results for the ensemble averaged population size are 
compared with simulations for single realizations of the process and also with standard mean field 
theory which is expected to be valid on large fully-connected graphs. 
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I. INTRODUCTION 



Infections produce further infections. This observation has long inspired theoreticians to 
find simple tractable evolution equations to model such a situation. One traditional and 
rather simple approach is the so-called SIR model originally introduced in 1], see also ^ 
and [3]. Here a certain population is divided into three distinct classes: the susceptible, S, 
wherein the individual is healthy but is allowed to catch the disease; further there are the 
infectious, denoted as I which is infected and can transmit the disease and the recovered R 
which is immune to further infection 4j 

Although the model is quite simple, it captures important features of the temporal dy- 
namics of an infection cycle. In so far the model is appropriate to describe a well-localized 
disease outburst. Due to the coupling of the three different groups S, I and R the process 
is non-linear. Furthermore, as long as there is an infectious population the system is in an 
nonequilibrium state not characterized by any physical a priory principle such as detailed 
balance. Despite its simplicity the SIR-model has not been solved exactly if fluctuations, 
which inevitably occur in a real system, are included in its description. In this paper we 
present such an exact solution using a mapping of the underlying Master equation into a 
quantum formulation. There appears a whole hierarchy of evolution equations for certain 
expectation values which can be closed and from which among other things the exact time 
evolution of the expected population size for each class can be extracted analytically in 
closed form. 

Our effort can be grouped in the permanent attraction exerted by modern biology and 
social science to understand the evolution of cooperative behavior. It is well known that 
in unstructured populations, natural selection favors defectors over cooperators. For that 
problem we also need the insight provided by mathematical tools. The SIR model offers a 



simple approach by a set of evolution equations 



H, S, H, E 



.Il0l| . To discuss the spreading 



of epidemics the SIR model can be implemented on a network 11,], which is further discussed 



in 
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3,Q. 



15] • The general scheme and the properties of networks are elaborated in detail 



16(| . In such a network approach, individuals are represented by nodes which are in either 
of the three states S, I, R. Contact between individual is modelled by links between the 
nodes. For maximal connectivity, where each individual is in contact with every other, i.e. 
for the fully connected graph, one expects the deterministic standard mean field equations 
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for the SIR model to be valid for large population size even if the infection and recovery 
are stochastic. In contrast, fluctuations are expected to cause substantial deviations from 
the mean field behaviour if the connectivity is low. Hence it is highly desirable to study 
the opposite case of minimal connectivity as realized in a linear chain. In the present paper 



we analyze the SIR model on a linear chain based upon the master equation 18] which is 
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22 1 . The method enables 



reformulated in terms of second quantized operators 
us to apply the algebraic properties of spin operators to determine a closed set of evolution 
equations for higher-order cluster functions. These cluster functions describe groups of 
susceptibles which can be infected from the boundary of the region. The time evolution for 
the cluster functions can be closed which makes the problem exactly solvable. This allows 
for a quantitative comparison with the mean-field solution as given by the deterministic 
standard SIR model and also with the random behavior of single realizations of the process 
obtained from Monte-Carlo simulation. 

The paper is organized as follows. In Sec. 2 we first define the standard deterministic 
SIR model and then introduce the stochastic dynamics that we consider to account for 
fluctuations. In Sec. 3 we describe the mathematical apparatus required for obtaining the 
exact results. This section can be skipped by readers not interested in the mathematical 
details. For an introduction into the quantum approach used there we refer to the reviews 



20, 



2l|. In Sec. 4 we present the exact results for the expected population densities and 
compare them with analytical results from the mean-field approach. In Sec. 5 we discuss 
results of Monte-Carlo simulations for single realizations of the process . In Sec. 6 we finish 
with some conclusions. 

II. STOCHASTIC SIR MODEL ON A GRAPH 

Let us denote with S(t), I(t) and R{t) the number of susceptibles, infectious and recovered 
individuals. The total number iV is conserved 

S(t) + I(t) + R(t) = N. (1) 
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In the traditional treatment of the SIR model the population strength is treated as a real 
number and infection and recovery are governed by the nonlinear set of coupled equations 



The first equation describes the decrease of the susceptible population through the infection 
of a susceptible individual by an infectious one. The loss is proportional to the infection rate 
(3 and since by definition S and I are non-negative, the loss is monotone. The second equation 
describes the gain of the infectious population of individuals by infection of susceptibles as 
described in the first equation and the spontaneous recovery with rate 7. The last equation 
follows simply from the conservation of the total number of individuals. 

These equations describe a deterministic evolution for each population class which entirely 
neglects fluctuations and which offers no description of the state of an individual member 
of the entire population. These equations may be regarded as a mean-field treatment of 
some fluctuating random process and therefore we shall refer to this standard SIR model 
as mean field SIR model. In view of our further approach it is appropriate to introduce 
the population densities nx(t) = X(t)/N where X stands for one of three classes S,I,R. 
Obviously the densities satisfy ns(t) + n/(t) + nn(t) = 1. 

We now define a stochastic SIR dynamics that describes the state of each individual. This 
description allows for randomness and hence fluctuations in both the infection and recovery 
process. In our individual-based version of the model the individuals are represented by the 
nodes of a graph. For each node % we introduce state variables which specify the state of the 
node. For reasons that become clear below we represent these state variables as occupation 
numbers which take value or 1 as follows: If node i is in the susceptible state at time t, 
we say that ns(i,t) = 1. If node i is in the infectious or recovered state then ns(i,t) = 0. 
Likewise we define occupation numbers nr(i,t),nR(i,t) which by definition are subject to 
the constraint ns(i,t) + rij{i,t) + nn(i,t) = 1. With this definition we define the (random) 
population sizes of class of individuals 
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where the sum is taken over all nodes of the graph. Considering N nodes ensures a strict 
conservation law S + I + R = N analogous to ([I]) for the deterministic SIR model. Contact 
between two individuals is represented by a link between two nodes. This defines the graph. 

The stochastic dynamics of the model is realized by the following Markov process. A 
susceptible individual at a node i becomes infectious after an exponentially distributed 
random time with rate /3I(i,t) where 

I(M) = X>&*) (4) 

iW 

is the total number of infectious individuals which are in contact with i at time t. This 
quantity is an integer random variable that depends on the current state of the system. On 
the other hand, an infectious individual at node i recovers after an exponentially distributed 
random time with fixed rate 7. Once an individual is recovered it remains so. All infection 
and recovery processes occur independently of each other. 

Thus this stochastic process is in double contrast to the evolution studied in the mean- 
field approach. There recovery and infection are deterministic and the infection rate is 
proportional to the size of the full population of infectious individuals. The latter property 
is recovered in our individual-based approach if each individual is in contact with every 
other, i.e., if the underlying network is the complete graph of N nodes. If then in addition 
the population size N is send to infinity, one expects fluctuations to disappear by the law 
of large numbers. Hence the traditional SIR model as described by (j2J) may be regarded as 
a deterministic limit of the evolution of our stochastic process on a complete graph in the 
thermodynamic limit of infinite population size. 

In our stochastic model the main quantity of interest is the expected state (nx(i,t) ) of 
a node % at time t, given some initial distribution. We shall focus on uncorrelated random 
initial distributions with some given mean population size for infectious and susceptible 
individuals. In this case the expectation value is independent of the node % and we write it 



in slight abuse of notation nxif) = (nx{i,t) ) and X(t) = (X(t) ) 23|. Moreover, in order 
to quantify and highlight the possible effect of fluctuations due to incomplete connectivity 
between individuals we study the most "non-mean-field" setting possible. I.e. we consider 
the lowest possible connectivity between individuals which is realized in a periodic chain 
of N nodes. The dynamics of the model is then realized by the following transitions on 
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neighboring nodes 



I S ^ 1 1 with rate (3 
S I 1 1 with rate j3 
I — > R with rate 7 . (5) 

The first two processes mean that a single susceptible S can catch the disease when it is 
in contact with infectious individuals situated on the neighboring nodes. The last process 
in Eq. ([5]) characterizes the recovering process, where an infectious individual recovers and 
becomes immunized, independently of the state of other individuals. 



III. QUANTUM APPROACH TO NONEQUILIBRIUM SYSTEMS 
A. Master equation in a quantum Hamiltonian representation 

Since the combined influence of noise and spatial degrees of freedom is an important 



issue in a theoretical understanding of biological and ecological processes [17| the dynamics 
due to Eq. ([5]) is formulated in a master equation for the full probability distribution of the 
process. Here we use a very transparent method, the transformation of the master equation 
into a quantum language. This exact mapping enables us to get an exact solution for the 
process defined above. Since the method borrows techniques from condensed matter and 
particle physics, we use well-established jargon that is slightly different from that above. In 
particular, we shall refer to nodes of the graph as sites on a lattice, and to the state variable 
nx{i) as occupation numbers by particles of type X. 

Let us summarize briefly the most important steps, for a detailed account of the approach 

sec 
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2lj . The master equation for the joint probability P(n,t) reads 

d t P(n,t) = CP(n,t) ■ (6) 

Here n stands for a certain configuration of S, I and R particles at time t. In a lattice gas 
description each lattice point is either empty or single occupied leading to nx{i) = 0, 1 for 
each type. Using the expansion 

\F(t)) = J2P(n,t)\n). (7) 



Eq. ([6]) can be rewritten as an equivalent equation in a Fock-space 

d t | F(t)) = —H | F(t)), (8) 

where the operator H is determined in such a manner that its matrix elements correspond 
to those of C. The formal solution of that equation is 

\F(t))=e- m \F(0)). (9) 

This expression gives the probability distribution at time t in terms of the initial distribution 
at time t — 0. 

It should be emphasized that the procedure is up to now independent of the realization 
of the basic vectors | n) . As shown by Doi [19] the average of an arbitrary physical quantity 
lZ{n) can be calculated by the average of the corresponding diagonal operator R{t) 

(R(t)) = p (n, t)K(n) = (s\R\ F(t)) (10) 

with the summation vector (s \= ^2(n |. The evolution equation for an operator R(t) can 
be cast in a commutator relation which reads 

d t (R) = (s\[R(t),H)_\F{t)). (11) 

As the result of the procedure, all the dynamical equations governing the classical problem 
are determined by the structure of the evolution operator H and the commutation rules of 
the operators. 

The evolution operator for the process defined by Eq. ([5]) reads 



- H = p J2 [ & l+i a *+i + h \-i a i-i ~ 4+i (1 - B i+1 ) - - Bi-i) 



B, 



+ 7 $>;-£,). (12) 

% 

Here a,, a\ and bi, b\ are the annihilation and creation operators for S and / types. The op- 
erators Ai = a\a>i and Bi = b\b{ represent the particle number operators with the eigenvalues 
and 1. The particle number operator A% corresponds to the occupation variable ns(i) and 
Bi corresponds to ni(i). 

The meaning of the evolution operator defined in Eq. (fl2|) is now transparent: The first 
term on the right hand side describes the annihilation of a susceptible at site i + 1 and a 



simultaneous creation of an infectious at the same site provided the neighboring lattice site 
i is occupied by an infectious indicated by the number operator Bi = b\bi. Mathematically 
this property is manifest in the commutator relation 

Mj] = (i - 2&t& i )<y y . (13) 

and similar rules for the a and a 1 . The operators commute at different lattice sites and 
the anticommute at the same lattice site. The anticommutator rule implies the exclusion 
principle, i.e. the eigenvalues of the particle operators are restricted to 0, 1 and therefore 
the corresponding averages fulfills < (Aj) < 1. Similar relations hold for Bi. From the 
definition follows 

(n s (i,t)) = (A) (14) 
(m(i,t)) = (B^ (15) 

and correspondingly 

(n R (i,t)) = l-(A i )-(B i ) (16) 

for the probability of finding node % in the recovered state. Notice that these expectation val- 
ues imply a double average over the initial distribution and over realizations of the stochastic 
dynamics. 

B. Cluster functions 

Using Eq. (PJ and Eq. (EDI we get 



= -!3{(A r B r _ l ) + (A r B T+1 )) 

at 

^(B r ) = P{(A r B r ^ l ) + (A r B r+l ))- 1 (B r ) (17) 

These equations involve second-order correlators which hints at the non-linear nature 
of the problem. To analyze the situation let us further study the higher order correlators 
appearing in Eqs. (TTTT) . For illustration we present the result for the two-point correlator 

^{A r B r+1 ) = -( 7 + P){A r B r+1 ) + /3 ((A r A r+1 B r+2 ) - (B r ^A r B r+1 )) (18) 



S 



To make a more systematic approach let us define the n-point cluster functions 

H r {n) = (A r A r+ i . . . A r+n _iB r+n ) 

Grin) = {B r _iA r A r+ i . . . A J . +n _iS r+n ) . (19) 

Obviously these functions are zero if one of the sites inside of the cluster is recovered with 
probability 1. Furthermore the functions H[n) and G(n) are sensitive to the fact that a 
cluster of susceptible individuals is diminished by infection at the border of the cluster. The 
introduction of these cluster functions is the decisive trick of our treatment which makes the 
nonlinear problem solvable. 

The cluster equations simplify under the natural assumption of a translation invariant 
initial distribution. Then one can drop the r-dependence and after a straightforward calcula- 
tion only taking into account the algebraic properties Eq. (fl3l) we end up with the following 
set of coupled equations for the cluster functions 

^H(n) = -( 7 + p)H(n)+/3(H(n+l)-G(n)) 

J^G(n) = -2(-f + (3)G(n) + 2f3G(n + l) for n > 1. (20) 

For non-translation invariant distributions the evolution equations for the cluster functions 
still close, but retain an extra r-dependence. The meaning of the evolution equations is 
immediately visible. The cluster described by H(n) decreases by infecting node r with rate 
(3 and recovery of node r + n with rate 7. Furthermore the H(n) cluster grows by increasing 
its length from n to n + 1 through infection of node n + r. The exact evolution equation of 
the G-cluster can be made plausible in similar terms. This competition between growth and 
reduction processes of the clusters yields a non-trivial steady state of the process discussed 
below. 

The second cluster equation can be solved recursively by treating the term G(n + 1) 
as an inhomogeneity of the remaining homogeneous first-order linear ordinary differential 
equation. One obtains 

G(n, t) = e~ 2 ^ + n > 0) . (21) 

1=0 L 

where G(n + I, 0) is an arbitrary initial condition. 
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Inserting that solution in the first equation we find in a similar fashion the solution 

H(n,t) = e-h + Mf^^H(!+n,0) 
1=0 l 

- P f dtfe-MW ^^ Gin + I, t') (22) 
^° 1=0 

These results are exact for arbitrary translation invariant initial distributions. For un- 
corrected random initial distribution of each class of individuals the initial conditions for 
the cluster functions read 

G(n, 0) = ns(0) n n/(0) 2 , H(n, 0) = n 5 (0) n n 7 (0). (23) 

Here the initial densities of infectious and susceptible individuals are ri/(0) and ns-(O). In- 
serting this in Eq. pi]) and Eq. ff22l results in 

2t 

G(n,t) = n 5 (0) n n 7 (0) 2 exp( ) (24) 

r 

H(n, t) = n s (0) n nj(0) exp(-- ) [1 - /3n r (0)r(l - exp(-t/r)] , (25) 

r 

where the relaxation time r is defined by 

T = 7 + /?(l- ns (0))- (26) 

We draw attention to the fact that the relaxation time depends on the initial conditions 
through the initial density of susceptibles. This is a consequence of the highly non-ergodic 
and far-from-equilibrium nature of the process. Both cluster functions decay monotonically 
in time to zero. 



IV. EXACT SOLUTION 

Using the exact result for the cluster expectation values for random initial conditions 
obtained in the previous section from the quantum approach, we are now in a position to 
obtain the exact time evolution for the expected number of individuals of each class. In 
terms of the cluster functions G(n) and H(n) the evolution equations for the quantities 

n s (t) = (A r ) and m{t) = (B r ) (27) 
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read 



^n s (t) = -2/3H(l,t) 
d 



dt 



m(t) = 2/3ff(l,t)- 7 n I (t). 



(28) 



Notice that H(l, t) is strictly positive for all finite times. Hence the density of susceptible in- 
dividuals is strictly monotonically decreasing which follows from the fact that no susceptibles 
are generated in the process. 

The averaged number of susceptible person at time t follows immediately from Eqs. (|24|) , 
leading to 



n s (t) = n s {0) - 2(3rn s {0)n I {0) 
+ 



1 < 2t \ 
1 — exp( J 

T 



1-exp — 

T 



[l-/3rn 7 (0)] 



(29) 



The decay of the susceptible persons S is not purely Arrhenius-like but given by a superpo- 
sition of two relaxation times. 

In the long time limit we find a nonzero stationary solution 



= n s (t -> oo) = n s (0) [1 - /3rnj(0)f 



(30) 



which can be written in the more transparent form 



715(0) 



j + Pil-nsW-miO) 



(31) 



7 + ^(1-715(0)) 

which makes the dependence on initial conditions and the recovery/infection ratio / j//3 fully 
explicit. 

In the same manner we find the expected density of infectious persons 



m(t) = nj(0)[exp(-7f) + 2prn s (0)f(t)} 



(32) 



with 



= eM - t /r)-eM- lt ) + MOlr ^ 

7r — 1 7r — 2 

Due to Eq. ([T]) 1 — n s — rij is the expected density of recovered individuals. 

For a comparison with the predictions of the original SIR model defined by the set of 
differential equations (j2J) the mean-field solution we highlight some features of this mean 
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field model. From the first equation in Eq. (j2j) we conclude that S(t) is a monotonically 
decreasing function. Moreover, the last equation implies that the stationary value for the 
infectious class is I* = 0. Both properties are shared by our stochastic SIR model. In 
the condition for the existence of a maximum in the number of infectious individual, the 
situation is more subtle. Writing the second equation in Eq. (J2D in the form 



one realizes that a maximum occurs if 5(0) > j//3. It is reached at a time where 
S(Q) = 7/ '/3. Inserting this into the first equation one may write in terms of the normalized 
population densities h s (Q) = 772/(6). Interestingly, the exact relation (128]) asserts that in 
our stochastic model the maximum occurs at a time determined by the same relation in 
the case of random initial conditions. Hence, our model reproduces several key features 
of the original SIR model. For these features, the low connectivity and stochasticity are 
unimportant for a comparison of average behavior of the stochastic dynamics with the 
deterministic behavior of the mean field model. Notice though, that the actual value of 
is not the same in the two models. It is also interesting to observe that the mean field 
expression for can be written in terms of the concentration in the form ng(G) = j/(/3N), 
i.e. the maximum in / occurs at a time where an initial concentration of susceptibles of 
order 1 has almost disappeared and only a finite total number (of order 1) of susceptibles 
are left in the population. This is in contrast to our stochastic model where the maximum 
in I occurs at a concentration of susceptibles which is of order 1. In this respect, the mean 
field model fails to capture the effects of low connectivity. 

For a more detailed analysis of the mean field SIR model we introduce g(t) = In S(t) 
which due to the first equation in Eq. ([2]) satisfies g = —(51. Differentiating again and using 
the second equation gives an integrable second order equation for g. After one integration 
one obtains 



where A is an integration constant. In case of the initial conditions S(0) = So, 1(0) = Iq 
and R(0) = it results A = 7 In So — f3(S + -M- Combining the last relation for g(t) with 
Eq. (|2D we then find the relation 



/ = IifiS - 7) 



(34) 



-f = (3e 9 - 7# + A 



(35) 



I(t) 



(3S(t) + j- In 




+ N. 



(36) 



12 



In the same manner we find 



m = -j^, (37) 



where R(t) obeys 



^ = - 7 <? exp[-^i?] + 7 (iV - R). (38) 
at 7 



It corresponds to an overdamped motion in a potential 

dR dU(R) . , rr/nx S 7 2 r m lR fr . n ^ 

This equation of motion does not allow for a closed solution in terms of elementary functions. 

In the limit t — > oo Eq. (I36I) gives a transcendental equation for the stationary population 
of susceptibles 

(3S* — In -— ^ — (5N = 0. (40) 
This has no solution in closed form, but for large N one obtains 

S* ~ S exp(--N) (41) 

7 

which decays exponentially in the population size N. This result is strongly different from 
the exact solution ( 13T|) where one finds a finite stationary value of order 1 even for infinite 
N. 

As a final remark we point out that in the mean field approximation one decomposes 
higher order correlators according to (AB) = (A)(B). Identifying (A) with the density of 
susceptibles n,s(t) and correspondingly (B) = ni(t) we get from (|28|) mean field equations of 
the form (jSJ), but with an infection rate f3 m f = (3/N. Hence the mean field approximation of 
our stochastic model yields a deterministic SIR dynamics with renormalized infection rate 

Pmf- 



V. MONTE CARLO SIMULATION DATA 



Our exact results are obtained for the thermodynamic limit of infinite population size, 
and they are results for a statistical ensemble of processes, averaged both over random 
initial states and histories. Here we present Monte Carlo simulation results for single runs 
of the process which demonstrate that even if the moderate population size is moderate, 
fluctuations around the computed expectation value are rather small. This mean that the 
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computed expectation values represent the typical behavior that one expects in a single 
outbreak of the disease. Only for very small populations the fluctuations around the expected 
mean become significant. 

We have performed the numerical simulation of the problem as follows. Initially, each 
site is occupied independently and randomly by a susceptible with probability n s (0) and 
by an infectious individual with probability n/(0) = 1 — n$(Q). For the dynamics we have 
chosen a random sequential update algorithm as follows. An arbitrary lattice site j is chosen 
randomly. If this site is occupied with an infectious I, then the /-state decays to R with 
probability 7/(2/? + 7). If it does not decay, then with equal probability 1/2 an adjacent 
site on the left or right hand site is chosen. If the chosen neighboring site is occupied by 
a susceptible S, then S is converted into I with probability 2/3(2/3 + 7). If lattice site j 
is occupied by a susceptible S or the site is recovered nothing happens. Then a new site 
is selected randomly and the procedure is repeated. N such update steps then define one 
Monte-Carlo time step. We remark that for an efficient implementation of the process one 
may keep a list of coordinates of infectious sites and select sites only from this list. However, 
for population sizes of the order of 10 3 such optimization is irrelevant for the study of single 
realizations of the process. 

In Fig. [1] we show simulation data for two different runs with population size 1000, 
demonstrating the absence (Fig. QJa)) or presence (Fig. [D(b)) respectively of a maximum 
in the number of infectious particles. The maximum occurs at values of n s which is of 
order 1, rather than 1/N as predicted by mean field theory. The finite limiting value of 
the susceptible population density is also clearly seen. The corresponding mean field value 
« 400 exp(— 7000) would be nearly zero. For a comparison of this single run with the 
computed mean values the corresponding exact expressions ns(t) (129]) and n.j(i) are shown 
as well. The deviations are at most in the range of a few percent. (132|) . 

In Fig. [2] we show simulation data and the exact solution for different total number of 
individuals. Fig. [2] demonstrates the he increasing effect of fluctuations for small population 
sizes. For population sizes of the order of 10 4 fluctuations become irrelevant. 
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5000 10000 15000 20000 40000 60000 80000 100000 120000 

t t 



(a) (3 = 0.7, 7 = 0.8, n s (0) = 0.4, nj(0) = 0.6 (b) (3 = 0.9, 7 = 0.1, n s (0) = 0.8, n,(0) = 0.2 

FIG. 1: Time evolution of susceptibles and infectious individuals for different rates (3 and 7. The full 
line corresponds to the simulations, the dashed line represents the exact solution. The population 
size is in both cases N = 10 3 




500000 1000000 1500000 2000 4000 6000 8000 10000 12000 14000 

t t 



(a) N = 10000 (b) N = 100 

FIG. 2: Time evolution of susceptibles and infectious individuals for different population size: The 
rates are in both cases (3 = 0.9, 7 = 0.1, the inital values are 77-5(0) = 0.8, nr(0) = 0.2. The dashed 
line is the exact solution. 
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VI. CONCLUSIONS 



In this work we have analyzed a SIR model for a population of susceptible S, infectious / 
and recovered R individuals evolving under a stochastic dynamics. In order to study the 
effect of fluctuations due to incomplete contact between the individuals we have defined the 
model on a linear chain. As an appropriate tool we have considered the master equation 
for the probability density which we wrote in a quantum formulation in terms of second 
quantized operators. These operators obey the commutation rules of Pauli operators, i.e., 
they commute at different lattice sites and anticommute at the same site. This property led 
us to find a coupled set of evolution equations forcertain cluster functions. These clusters 
describe the behavior of susceptibles surrounded by infectious individuals at the edges of 
the clusters and allow for an exact analytical treatment of the whole hierarchy of evolution 
equations. We stress that in the exact solution all fluctuations are included. 

Comparing this exact solution with the behavior of the traditional mean field SIR model, 
we find a significant difference. Whereas the mean field solution yields a stationary density 
for the susceptibles n* s ~ exp(— f3N/^) which depends on the population size and is extremely 
small for large N, the exact solution reveals a stationary density independent of N and of 
order 1. This shows on a quantitative level how fluctuations due to low connectivity of 
individuals are crucial for understanding the spreading of a disease in the framework of the 
SIR mechanism. 

We remark that by making a mean field approximation to the exact evolution equations 
(I28p of our model, one obtains a deterministic set of evolution equations similar to those 
of the mean- field SIR model, but with an infection rate j3 m f = j3/N . Indeed, inserting /3 m f 
in the stationary density of the mean field SIR model, yields a finite stationary density of 
susceptibles of order 1, as in our stochastic SIR model. Thus the effect of low-connectivity 
model can be qualitatively described by a mean-field model with a small renormalized in- 
fection rate j3 m f. Capturing the precise form of the time evolution, however, is beyond the 
capabilities of the mean-field description. 

The analytical findings are illustrated by numerical simulations which demonstrate that 
fluctuations due to finite population size are negligible for population sizes of order 1000 or 
larger. We stress that while here we have focused on uncorrelated initial distributions which 
are on average spatially homogeneous, our exact analytical approach can be extended to 
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study the effect of correlations and spatial inhomogeneities in the initial distribution. The 
model remains exactly solvable also for finite population size. 
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